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Abstract. We consider the extended Hubbard model in the atomic limit on a Bethe lattice with coordina- 
tion number z. By using the equations of motion formalism, the model is exactly solved for both attractive 
and repulsive intersite potential V. By focusing on the case of negative V, i.e., attractive intersite inter- 
action, we study the phase diagram at finite temperature and find, for various values of the filling and of 
the on-site coupling U, a phase transition towards a state with phase separation. We determine the critical 
temperature as a function of the relevant parameters, i//|V|, n and z and we find a reentrant behavior in 
the plane (U/\V\,T). Finally, several thermodynamic properties are investigated near criticality. 

PACS. 71.10.Fd Lattice fermion models - 71.10-w Theories and models of many-electron systems 



1 Introduction 

Understanding the effects of competing interactions and 
the corresponding quantum phase transitions in models of 
strongly correlated electron systems is a crucial problem 
in condensed matter physics. A thorough comprehension 
could shed new light on the behavior of new classes of 
novel materials, ranging from inorganic chain compounds 
and conducting polymers to organic and high temperature 
superconductors. In all such systems, long-ranged Cou- 
lomb interactions play an essential role, as it has been 
pointed out in Ref . [1] , where it was shown that the simple 
Hubbard model does not offer a model for the supercon- 
ductivity in the right order of amplitude. Similarly, Frolich 
and Coulomb long-ranged interactions have been recog- 
nized to play an essential role in a realistic multi-polaron 
model of high-T c superconductivity [5] . 

Here, we study the so-called extended Hubbard model 
(EHM) [3] , where a nearest-neighbor interaction V is added 
to the original Hubbard Hamiltonian, containing only an 
on-site interaction U. The intersite interaction can some- 
how mime longer-ranged Coulomb interactions. In various 
applications of the model, the parameters U and V can 
represent the effective interaction couplings taking into ac- 
count also other interactions (for instance with phonons). 
Therefore we assume that U and V could take positive as 
well as negative values. The EHM has been deeply investi- 
gated by means of various analytical techniques. In partic- 
ular, g-ology investigations [I] as well as renormalization 
group [5] and bosonization 6 analysis have been success- 



fully carried out in the weak coupling regime, while exact 
diagonalization calculations [7 and quantum Monte Carlo 
simulations [8] have been used to study the intermediate 
and strong coupling regimes [9]. The cases of half filling 
(n = 1) 8J and quarter filling (n = 0.5) [TU] have been 
deeply investigated in order to establish the boundary 
line between spin-density wave and charge-density wave 
phases, as well as existence of an intermediate bond-order 
wave phase. More recently several systematic studies of 
the phase diagram of the one-dimensional extended Hub- 
bard model as a function of all the three parameters n, 
U, V have been performed as well [llj . revealing a very 
rich structure due to the presence of multiple first and sec- 
ond order phase transitions, critical points and interesting 
reentrant behaviors [12] . The relation with the experimen- 
tal findings relative to some manganese compounds [13] 
has been investigated in detail together with the possibil- 
ity of improving the description of cuprates with respect 
to the simple Hubbard model [12] . However, despite of 
the above results, such a general case still lacks a detailed 
study. 

In order to overcome the difficulties shown by the EHM 
and to investigate charge ordering in interacting electron 
systems, one may consider a minimal model, allowing one 
to capture all the relevant phenomenology while being ex- 
actly solved. That is, one can consider the very narrow 
band limit, where tij <C U, V; for simplicity the classical 
limit = can be taken. In this way one is led to the 
extended Hubbard model in the atomic limit (AL-EHM), 
which is one of the simplest models apt to describe charge 
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Fig. 1. Distribution of the particles along the Bethe lattice at 
n = 0.5 and T = 0: (a) V < and U > 2V; (b) V < and U < 
2V. White, grey and black squares denote empty, arbitrary 
spin singly occupied and double occupied sites, respectively. 



ordering in interacting electron systems. The charge or- 
dering has attracted much attention after the discovery 
of charge-density waves in high temperature superconduc- 
tors and since then the problem of competition and coexis- 
tence of this kind of ordering and other orderings has been 
widely investigated [14] . Charge ordering is in fact accom- 
panied by metal-insulator transitions through the growth 
of charge-order parameter and the opening of gaps at the 
Fermi surface. Furthermore, phenomena such as colossal 
magnetoresistance in perovskite-type manganese oxides 
are observed immediately when the charge order melts [T5] 
ITS] . Recently charge orderings and charge-density waves 
have been widely observed in transition metal compounds 
[17] and also in organic conductors [18]. It is worth notic- 
ing that the observed changes in transport, optical, di- 
electric and magnetic properties at or near charge-order 
transitions have spurred a lot of interest on the problem of 
their control through the charge-order transitions. Thus, 
in this context, a deeper understanding of the nature of 
charge-order transitions could be relevant as well as of the 
whole structure of the phase diagrams. Another aspect 
which is very interesting and deserves further investiga- 
tions is the phase separation phenomenon p!9 ll20ll2T] : it 
has been investigated in some detail in the AL-EHM con- 
text by means of Monte Carlo simulations and Hartree- 
Fock approximation in Ref. [22] and the results have been 
compared with the experimental ones obtained for the or- 
ganic conductor (DI — DCNQI)2Ag [23]. The phase sep- 
aration phenomenon is characterized by a ground state 
macroscopically inhomogeneous with different spatial re- 
gions with different average charge densities. Typically, 
it occurs for attractive on-site and intersite interactions, 
so that the electrons cluster together in order to gain the 
maximum negative potential energy. In the extended Hub- 
bard model, two types of phase separated configurations 
are possible: a cluster of singly occupied sites or a clus- 
ter of doubly occupied sites, as shown in Fig. [T] A cluster 
of doubly occupied sites is favored when the on-site and 
intersite interactions are both attractive, whereas cartoon 
(b) would correspond to large positive U and negative V. 

The critical behavior of the AL-EHM has been stud- 
ied by many authors, mainly in the one dimensional case 
[2^[2^[2~rJl[27P51[2^[3^[3T|l3iT^ . and a staggered charge 
order has been found irrespectively of the method of anal- 



ysis. Very recently, a comprehensive and systematic anal- 
ysis of the one dimensional AL-EHM has been carried out 
both at T = and finite temperature in the whole pa- 
rameters space and various types of long-range charge or- 
dered states have been observed [33] . A study of the two- 
dimensional version of the model has been performed [20] 
[35] . and very recently the same model on a Bethe lattice 
with general coordination number z has been investigated 
[3J. 

The study of the extended Hubbard model with attrac- 
tive as well as repulsive on-site and intersite interactions 
could also shed new light on the properties of systems of 
ultracold fermionic atoms loaded in optical lattices, such 
as fermionic superfluid properties in a two dimensional 
lattice [37] , and the possibility of the supersolid state [38] . 

In this paper we study the AL-EHM on a Bethe lat- 
tice with coordination number z within the equations of 
motion framework [33 l l39 l l40] . This formalism allows us to 
exactly solve the model. By focusing on the case of attrac- 
tive V, we study the phase diagram at finite temperature 
T for different values of the filling and the coupling U/\V\. 
By varying the temperature, we find a region of negative 
compressibility, hinting at a transition from a thermody- 
namically stable to an unstable phase characterized by 
phase separation. The critical temperature, function of the 
external parameters U/ \ V\, n and z, shows a reentrant 
behavior in the plane (U / \V\ , T) at half filling (n = 1). 
Although thermodynamically unstable, the study of the 
T < T c phase allows us to determine the critical value 
of the on-site interaction U c separating the formation of 
clusters of singly and doubly occupied sites. 

The plan of the paper is as follows: in Sec. [2] we briefly 
review the equations of motion method for the extended 
Hubbard model in the atomic limit on a Bethe lattice with 
coordination number z and compute the Green's and cor- 
relation functions. We find that they crucially depend only 
on two parameters, which can be self-consistcntly com- 
puted allowing us to determine all the local properties of 
the system in the case of general filling n. In Sec. [3] we de- 
rive the phase diagram in the planes (n, T) and (U, T) for 
a Bethe lattice with coordination number z = 3, and we 
show the behavior of the relevant thermodynamic quanti- 
ties. Finally, Sec.[4jis devoted to our concluding remarks. 

2 The general model 

The extended Hubbard Hamiltonian in the narrow-band 
limit on a Bethe lattice with coordination number z and 
nearest-neighbor interactions is given by: 

z 

H = -H0) + C/D(0) + ^i2 (p) . (1) 
P =i 

Here represents the Hamiltonian of the p-th sub-tree 
rooted at the central site (0): 

z-l 

ff&O = -nn(p) + UD(p) + Vn(0)n(p) + ^ H^< m \ (2) 

m—1 
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where p (p = 1, . . . , z) are the nearest neighbors of the site 
(0), also termed the first shell. In turn, ij(P' m ) describes 
the TO-th sub-tree rooted at the site (p) , and so on to infin- 
ity. U and V are the strengths of the local and intersite in- 
teractions, /i is the chemical potential, n(i) = n|(i) + rtj(z) 
and D(i) — n^(i)ni(i) = n (i) [n (i) — 1] /2 are the charge 
density and double occupancy operators at a general site 
i . As usual, n a (i) = c\,(i)c a (i) with a = {|, j} and c a (i) 
(cj-(z)) is the fermionic annihilation (creation) operator of 
an electron of spin a at site i, satisfying canonical an- 
ticommutation relations. We adopt the Heisenberg pic- 
ture i = (i,t), and we use the spinorial notation for the 
fermionic fields. 

The exact solution of the model can be obtained by us- 
ing the equations of motion approach in the context of the 
composite operator method [41]. Following this line, the 
main point is the choice of the field operators. A conve- 
nient operatorial basis is given by the Hubbard operators, 
= [n(i) — l]c(i) and r/(i) = n(i)c(i), which satisfy the 
equations of motion: 



d 

l Ql r l(' 1 ) = \U ~ tivW + zVr)(i)n a (i). 



(3) 



Hereafter, for a generic operator 4> (i) we shall use the no- 
tation 4> a (i) = j^* =1 (hp) being the first near- 
est neighbors of the site i. The Heisenberg equations ((3|) 
contain the higher-order nonlocal operators £,(i)n a (i) and 
r](i)n a (i). By taking time derivatives of the latter, higher- 
order operators are generated. This process may be con- 
tinued and usually an infinite hierarchy of field operators 
is created. However, by noting that the number n(i) and 
the double occupancy D(i) operators satisfy the following 
algebra 



riP (i) — n(i) + a p D(i), 
D p (i) = D(i), 
n p (i)D{i) = 2D(i) + a p D(i), 



(4) 



where p > 1 and a p = 2 P — 2, it is easy to establish the 
following recursion rule: 



[n«(i)] k = £A«[n«( 



(5) 



where 



satisfy the equations of motion 



d 



( \ 

I r]{i)[n a {i)] > 



\»7(»)[n«(*)p/ 
(7) 



(8) 



Here and are the energy matrices of rank (2z + 
1) x (2z + 1) [36] whose eigenvalues are: 



E$ = -M + (m - 1)V, 
E%>=-n+U+(m-l)V, 



(9) 



where m — 1, . . . , (2z + 1). The model has now been for- 
mally solved since one has found a closed set of eigen- 
operators and eigenvalues allowing one to ascertain exact 
expressions of the retarded Green's function (GF) 

G^(t - = 0(t - t')({^(0,t),r f (0,t')}), (10) 
and, consequently, of the correlation function (CF) 

C {s \t-t') = (^W(0,i)^W t (0,i / )>- (11) 

Indeed, by means of the equations of motion ([8]) one finds 
[36]: 



2z+l 



G«( W )=E 



(s,m) 



I U> — E„ 



iS 



(12) 



and 



CW(o;) = tt <T {s ' n) T£ S (w - E$) , (13) 

m— 1 

where s = £, r?, = 1 + tanh(/3E^} /2), (3 = l/k B T and 
Em are given in Eq. The spectral density matrices 
a ob can ^ e computed by means of the formula: 

2z+l 



r (s,n) _ 



an / ^ r; 
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(14) 



The coefficients Am are rational numbers, satisfying the 

relations J2m=i A $ = 1 and = <W (A = 1, • • ' . 2z ) 
[33] , The recursion relation (J5j) allows one to close the 
hierarchy of equations of motion. Thus, one may define a 
new composite field operator [36] : 



(6) 



In Eq. [[13]). ftM is the (2z + 1) x (2z + 1) matrix whose 
columns are the eigenvectors of the energy matrix 
and 1^ is the (2z + 1) x (2z + 1) normalization matrix 

whose elements can be written as [3H]: In,m = ft(™+ m ~ 2 ) — 
A („+m-2) j j(v) m = A („+m-2)_ Here the cnarge correlators 

O) and AW are defined as 

K W = ([n"(0)] fc ), A< fc ) = J(n(0)K(0)] fe ). (15) 
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By exploiting the recursion relation ((5]), it is not difficult to 
show that also and A^-* obey to similar recursion rela- 
tions, limiting their computation to the first 2z correlators 
[36] . However, the knowledge of the GFs and of the CFs, is 
not fully achieved because they depend on the external pa- 
rameters n, T, U, V, as well as on the internal parameters 



f, ...,2z). The internal parameters can 



fi, kW, \W (k 
be self-consistently computed by using algebra constraints 
and symmetry requirements [36]. Within this scheme, in 
Ref. [36] it has been shown that the charge correlators K^ k ' 
and A<*> can be written as a function of only two parame- 
ters, X\ and X 2 , in terms of which one may find a solution 
of the model. X\ and X2 are parameters of seminal im- 
portance since all correlators and fundamental properties 
of the system under study can be expressed in terms of 
them. Upon requiring translational invariance, one finds 
two equations allowing one to determine X\ and X 2 as a 
function of the chemical potential [36] : 

X 1 = 2e^(l - X 1 - dX 2 )(l + aXi + a 2 X 2 ) z ~ 1 



e m»-u) [2 + (d - \)X l - 2dX 2 ]{l + dX x + d 2 X 2 ) z ~ 1 

x 2 = e ^- u) [i + dx 1 - {2d + i)x 2 }{\ + dx x + d 2 x 2 y- 1 



2ePvK 2 X 2 (l + aXi + a 2 X 2 



(16) 



Here we defined K = ', a = (K-l) and d = (K 2 -l). 
In turn, the chemical potential fi can be determined by 
fixing the particle density: 

_ (X 1 -2X 2 )(l+aX 1 +a 2 X 2 )+2X 2 (l+dX 1 +d 2 X 2 ) 

11 ~ (l-X 1 +X 2 ) + (X 1 -2X 2 )(l+aX 1 +a*X 2 )+X 2 (l+dX 1 +d*X 2 ) ' 

(17) 

Equations fl6]) and (fT7]) constitute a system of coupled 
equations ascertaining the three parameters /1, X\ and 
X 2 in terms of the external parameters of the model n, U, 
V , T. Once these quantities are known, all local properties 
of the model can be computed. 



3 Results and discussion 

According to the sign of the intersite potential V, the solu- 
tions of the self-consistent equations (fTB]l and (fTT)) can be 
remarkably different. Here we study the case of attractive 
nearest neighbor interaction V (V < 0), allowing the on- 
site interaction to be both repulsive and attractive. Upon 
fixing V = — 1 and taking \ V\ = 1 as the unit of energy, we 
study the equations (fl6]l and (fT7]) for various values of the 
external parameters n, T and U, focusing on the simpler 
case z = 3. By considering higher coordination numbers, 
one obtains similar results. In Ref. [36] . upon fixing the 
chemical potential at the half filling value \i — zV + U/2, 
we have shown that a phase separation exists only for 
V < 0. Since experimentally it is possible to tune the den- 
sity in a controlled way, by varying the doping, here we 
fix the particle density n: the chemical potential will be 
determined by the system itself according to the values of 
the external parameters. 




Fig. 2. The chemical potential jj, as a function of the particle 
density n for V = -1, z = 3, T = 0.3, 0.7, 1, 1.5 and for (a) 
U = -2, (b) U = -0.5, (c) U = 0.5, and (d) U = 2. 



In Figs. [2^,-d we plot the chemical potential fj, as a 
function of the particle density n for both attractive and 
repulsive on-site interactions, and for several values of the 
temperature. In the high temperature regime, fi is always 
an increasing function of the particle density. By lower- 
ing the temperature, one observes that there is a critical 
value of the temperature below which the chemical po- 
tential becomes a decreasing function of n in the interval 
n-i < n < n 2 . The width of the interval depends on both T 
and U. In Figs. [3] we plot the derivative dn/dfi as a func- 



- T=0.3 
■■ T=0.7 
-T=1.0 

-T=1.5 





(a) 



(b) 





(c) 



(d) 



Fig. 3. The derivative dn/dfi as a function of the particle 
density n for V = -1, z = 3, T = 0.3, 0.7, 1, 1.5 and for (a) 
U = -2, (b) U = -0.5, (c) U = 0.5, (d) (7 = 2. 
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(c) 

Fig. 4. The critical temperature T c as a function of the par- 
ticle density n for V = —1, z = 3 and for (a) U = —2, -0.5, 
0.5, 2; (b) U — 2.3, 2.5, 3. (c) Comparison between the critical 
temperatures for z — 3 and z = 4. 



tion of n for the same values of T and U given in Figs. [H 
One observes that, at high temperatures, dn/d/j, is a reg- 
ular function, always positive, for all values of n. At low 
temperatures, dn/dfj, exhibits a divergence at the critical 
points n\ and 71.2, is negative in the interval ri\ < n < 712, 
and positive outside. Since the derivative dn/dfi is pro- 
portional to the thermal compressibility k = l/n 2 - dn/d/j,, 
one immediately infers that m < n < is a particle den- 
sity region where the system is thermodynamically unsta- 
ble. A possible solution to this problem is to abandon the 
translational invariancc assumption made in order to ob- 
tain the self-consistent equations (fT6|) . Nevertheless, here 
we are interested only in describing the transition to a 
phase separated regime, and, thus, we shall not address 
the problem of finding a thermodynamically stable solu- 
tion. By studying the temperature dependence of n\ and 
?i2, one can obtain the phase diagram in the plane (n, T). 
For a fixed value of n, there is a critical temperature T c 
at which the system undergoes a phase transition from a 
thermodynamically stable phase to an unstable one. Be- 
low T c the homogeneous state is not allowed due to a 
negative compressibility, and a phase separation occurs. 
In Figs. BJi-b we plot T c as a function of the particle den- 
sity for several values of U. The critical temperature T c 
and the width of the instability region An = — n\ in- 
crease by decreasing U: an attractive on-site potential will 
favor phase separation and in particular the clustering of 
doubly occupied sites. As a function of the particle den- 
sity, the critical temperature shows a lobe-like behavior: 
it increases by increasing n up to half filling, where it has 
a maximum; further augmenting n, it decreases vanishing 
at n = 2. A different behavior is observed when U > 2\V\ 
(see Fig. 0t>): the lobe splits in two parts and T c increases 
with n up to quarter filling, then decreases with a min- 



imum at half filling. The n = 1 critical temperature is 
finite only for 2V < U < U c , where U c = 2.3V for z = 3. 
For larger on-site potentials, there is no transition at half 
filling. For higher coordination number, one notices an in- 
crease of the critical temperature, as evidenced in Fig. 2J:, 
where we compare T c at U = — 2 for z = 3, 4. A useful 
representation of the phase diagram is obtained by plot- 
ting the phase boundary line in the plane (T c , U). In Fig. 
[5] we plot the critical temperature T c as a function of the 
on-site potential U for several values of n. Starting from 
large attractive U, T c decreases by increasing U and in- 
creases by increasing n, as already evident from Fig. [4] 
More interestingly, at half filling one observes a reentrant 
behavior with a turning point at U = U c w 2.3V (for 
z = 3) and at U — 2\V\ the critical temperature van- 
ishes. For particle densities less than half filling, one does 
not observe a reentrance in the phase diagram, but rather 
a flattening of T c at a finite value. Although below the 



2.0 




Fig. 5. The critical temperature T c as a function of the on- 
site potential U for V = — 1, z — 3 and for n — 0.25, 0.5, 0.75, 
1. 

critical temperature T c the system is thermodynamically 
unstable, the study of the behavior of relevant parameters 
below T c , unveils the existence of a critical value of the on- 
site potential separating the two possible realizations of a 
phase separated region, sketched in Fig Q] In Figs. [6] we 
plot the double occupancy D, the short-range correlation 
function A*- 1 ) and the internal energy E as functions of U 
at half filling and various temperatures. One observes that 
D and A^ 1 ^ exhibit two plateaus at low temperatures. In 
the limit T — ► there is a discontinuity around U c ~ 2.3V 
(for z — 3). For U > U c , the double occupancy tends 
to zero as T — > 0, whereas \^ tends to 1/2. The repul- 
sion between the electrons on the same site and the con- 
comitant nearest neighbor attraction, leads to a scenario 
where the electrons tend to cluster together singly occupy- 
ing neighboring sites. As a consequence, at half filling all 
sites are singly occupied, whereas for n < 1 one observes 
two separated clusters of singly and empty sites. When 
U < U c , one observes a dramatic increase (step-like) of the 
double occupancy and of the short-range correlation func- 
tion, namely: D — > 1/2 and A^ 1 ' — > 1. As a consequence, 
half of the sites are doubly occupied, and AW = 1 indi- 
cates that the electrons tend to occupy nearest neighbor 
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1 1.5 2 2.5 3 

U 



(c) 

Fig. 6. (a) The double occupancy D, (b) the short-range corre- 
lation function A' 1 ' and (c) the internal energy E as functions 
of U for V = -1, z = 3, n = 1 and for T = 0.15, 0.2, 0.3, 0.4, 
0.5. 

sites, arranging to form large domains occupied, leaving 
the rest of the lattice empty. For U < U c one observes 
also a dramatic decrease of the internal energy E, with a 
discontinuity as T — > around U ~ U c . In Fig. [7] we plot 
the specific heat C as a function of the temperature T at 
half filling and for several values of the on-site potential U. 
One observes that, for U — » U c , the peak becomes sharper 
and a divergence is observed at U = U c , confirming the 
emergence of a critical point. The same analysis carried 




(a) (b) 

Fig. 7. The specific heat C as a function of the temperature 
T for V = -1, z = 3, n = 1 and for (a) U = -1.0, -0.5, 0, 0.5, 
1.0; (b) (7 = 2, 2.1, 2.2, 2.25. 

out for different values of the particle density, evidences a 
similar behavior of D, A^ and E as a functions of U. 



4 Conclusions 

In this paper we have obtained the finite temperature 
phase diagram of a system of fermions with on-site and 
nearest neighbor interactions localized on the sites of the 
Bcthe lattice. The Hamiltonian describing such a system 
defines the so-called extended Hubbard model in the atomic 



limit. By means of the equations of motion method, it is 
possible to exactly solve the model. For attractive near- 
est neighbor interaction, by varying the temperature we 
found a region of negative compressibility, hinting at a 
transition from a thermodynamically stable to an unstable 
phase characterized by phase separation. The critical tem- 
perature T c at which the transition occurs depends on the 
filling n, on the ratio ?7/|V| and on z. It increases with in- 
creasing z and presents a reentrant behavior at half filling 
for J7/|V| > 2. Even if the region below T < T c is thermo- 
dynamically unstable, the study of several thermodynamic 
quantities allowed us to determine the critical value of the 
on-site interaction U c separating the formation of clusters 
of singly and doubly occupied sites. It would be interest- 
ing to investigate in detail the structure of the unstable 
phase as well as to look for a thermodynamically stable 
phase. This could be done by modifying, for instance, the 
translational invariance assumption made to obtain the 
self-consistent equations (fTo) and (fTT)) . 
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